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ABSTRACT 

A  family  of  iterative  schemes  is  presented  for  the 
solution  of  systems  of  elliptic  difference  equations.  The 
family  is  defined  by  a  generalization  of  the  usual  notion 
of  extrapolation  or  over-relaxation  and  includes  many  well 
knovm  procedures. 

Theorems  about  the  relationships  between  the  eigenvalues 
of  some  of  these  schemes  and  their  dependence  upon  the 
"extrapolation"  parameters  are  proved.  These  results  are 
used  to  define  three  subclasses,  called  complete  image 
classes,  of  the  general  family  of  schemes.  The  complete 
image  classes  are  such  that  each  eigenvalue  of  any  scheme 
in  a  class  is  a  given  algebraic  function  of  at  least  one 
eigenvalue  of  any  other  scheme  in  the  same  class.  Thus  a 
knowledge  of  the  eigenvalues  of  any  scheme  in  a  given  complete 
image  class  enables  the  determination  of  the  "best"  scheme  in 
that  class. 

Some  special  equations  are  considered  in  which  the 
eigenvalues  of  any  scheme  in  the  general  family  can  be 
obtained  in  terms  of  the  eigenvalues  of  two  matrices. 

All  of  the  results  are  applied  to  the  case  of  the 
Laplace  difference  equations,  as  an  example.  Some  well- 
known  results  on  convergence  of  familiar  schemes  are 
immediate  consequences.  Analogous  results  are  found  for 
the  less  familiar  "line"  or  block  methods. 
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NYO-7971 

ON  SOME  ITERATIVE  METHODS  FOR 
SOLVING  ELLIPTIC  DIFFERENCE  EQUATIONS 

1.  Introduction. 

The  nximerlcal  solution  of  boundary  value  problems 
for  partial  differential  equations  usually  requires  the 
solution  of  large  systems  of  linear  equations.  The  order, 
n,  of  such  systems  Is  essentially  equal  to  the  number  of 

mesh  points  In  the  domain  under  consideration.  Since 

3 
direct  Inversion  procedures  require  the  order  of  n"^ 

operations  they  are  not  practicable,  even  using  high  speed 

digital  computers,  for  reasonable  meshes  In  two  or  more 

dimensions.  Thus  Iterative  methods  for  solving  linear 

systems  are  of  great  Interest  as  they  usually  require 

the  order  of  n   operations.  In  addition  the  coefficient 

matrix  of  the  system  which  results  from  the  finite 

difference  approximations  has  many  strategically  placed 

zeroes.  However,  no  special  accoxont  of  these  zeroes  Is 

taken  In  most  direct  Inversions  or  In  general  Iterative 

procedures.  It  Is  reasonable  to  expect  that  particular 

methods,  designed  In  accordance  with  the  general  structure 

of  the  coefficient  matrix,  could  further  reduce  the 

number  of  operations.  Many  such  special  Iteration  schemes 
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have  been  devised  and  conditions  on  the  coefficient  matrix 
which  are  sufficient  to  Insure  the  convergence  of  some  of 
these  methods  have  been  obtained  [1,7,9].   However  there 
Is  no  general  comparison  procedure  to  determine  which 
of  many  possible  methods  Is  "best"  In  a  given  case. 

In  the  present  paper  we  formulate  a  family  of  Iterative 
schemes  for  a  particular  class  of  coefficient  matrices 
(in  which  the  zeroes  are  placed  as  In  the  usual  five- 
point  Laplace  difference  equations).   This  family  Is  de- 
fined by  a  generalization  of  the  usual  notion  of  extrapo- 
lation or  over-relaxation.   It  Is  then  possible  to  formulate 
the  problem  of  finding  the  "best"  scheme  and,  more  Important, 
some  general  theorems  on  the  eigenvalues  of  these  schemes 
are  proved. 

The  theorems  are  used  to  define  three  subclasses, 
called  complete  Image  classes,  of  the  general  family. 
These  classes  contain  many  of  the  schemes  In  current  use 
as  well  as  generalizations  of  them.   Thus  It  Is  shown 
that  a  variety  of  Independently  proposed  and  seemingly 
unrelated  Iterative  methods  are  special  cases  of  a  general 
class  of  methods.   These  complete  Image  classes  are  such 
that  each  of  the  eigenvalues  of  any  scheme  In  the  class 
Is  a  given  algebraic  function  of  at  least  one  of  the 
eigenvalues  of  a  particular  reference  scheme  of  the  class. 
Thus  a  knowledge  of  the  eigenvalues  of  the  reference  scheme 
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permits.  In  principle,  the  determination  of  the  best  scheme 
of  the  given  class. 

A  special  class  of  eqxiations  is  considered  for  which 
all  the  eigenvalues  of  each  reference  scheme  can  be 
explicitly  written  in  terms  of  the  eigenvalues  of  two 
matrices.   It  is  then  possible  to  determine  which  of  the 
reference  schemes  is  best  and  hence  it  is  also  possible 
to  determine  that  scheme  which  is  best  of  all  those  in 
the  complete  image  classes. 

As  an  example  of  the  application  of  this  general 
theory  the  Laplace  difference  equations  are  considered. 
The  well-known  results  on  rate  of  convergence  for  the 
Richardson,  Liebmann  and  extrapolated  Liebmann  methods 
are  immediate  consequences.   Analogous  results  are 
obtained  for  the  less  well-known  "line"  methods  which 
are  shown  to  be  superior. 

It  is  clear  that  the  methods  of  the  present  paper 


The  origin  of  iterative  line  methods  is  obscure.  They 
have  been  in  use  by  Russian  mathematicians  for  a  nijmber 
of  years.   About  19^5  J.  von  Neumann  and  L.  H.  Thomas 
independently  proposed  line  methods  for  parabolic 
difference  eq\aations.   An  independent  investigation  was 
initiated  by  M.  E.  Rose  and  the  present  author  in  1953 
and  some  of  the  results  of  Section  8  were  then  obtained. 
Peaceman  and  Ratchford  have  studied  their  applicability 
to  parabolic  difference  equations  and  double-sweep 
iterative  solutions  of  the  Laplace  difference  equations. 
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can  be  applied  to  more  general  Iterative  schemes  than 
those  considered.  In  addition  many  of  the  present  results 
can  be  easily  extended  to  partial  difference  eq\iatlons  In 
higher  dimensions  and  with  more  general  boundaries. 
2.   Formulation. 

A  large  class  of  two  dimensional  linear  elliptic 
difference  equations  are  of  the  form: 

(2.0)       -  'lA.l,J  -  ^/lJ-1  +  ^J  -  ^IJ^+IJ 

~  ^Ij'^lJ+l  "  ^IJ 

within  a  coordinate  rectangle  specified  by  1  <  1  <  p, 
1  <  J  <  q.   On  the  boundaries  of  the  domain  eqixatlons 
of  the  form  (2.0)  hold  with  the  coefficients: 


(2.1)     /^j  =  Tpj  «  0,  1  <  J  <  q; 

^11  "  *lq  '^  °'  1  <  i  <  P- 
Such  equations  are  obtained  from  second  order  elliptic 


The  notation  used  for  the  coefficients  Is  partly 
mnemonic.  In  that  the  letters  stand  for  "left", 
"bottom",  "right",  and  "top".  Indicating  which 
neighbor  of  the  central  mesh  point  Is  Involved. 


« 


The  five  point  scheme  Implied  by  (2.0)  assumes  no 
mixed  partial  derlvates  In  the  equations. 
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partial  differential  equations  by  applying  the  U8\ial 
second  order  difference  approximations  on  some  coordinate 
mesh,   (^»T|.  ).   The  coefficient  matrix  of  the  resulting 
system,  (2.0-1),  must  be  non-singular  and,  with  a  little 
care  In  differencing  [1],  can  be  made  positive  definite 
(and  symmetric  If  the  eq\aatlon  Is  self  adjoint).  However, 
we  shall  assiane,  unless  otherwise  stated,  only  the  non- 
slngularlty. 

For  convenience  of  notation  and  discussion  we 
Introduce,  for  each  J,  the  p-dlmenslonal  column  vectors: 


^ij 


(2.2) 


%- 


0. 


2J 


'J  - 


and  the  p  x  p  order  matrices: 


•*2J 


8 


PJ 


y 
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/ 


^ 


'2J° 


hi". 


•  • 


0  r 


^ 


'J  " 


V 


0  r 


2J 
0 


f      • 


,  1  <  J  <  q; 


.  ''p-iJ 


V 


(2.3) 


'  J 


^ 


'IJ       0 
^2J 


'j  - 


'PJ 


r 


,  1  <  J  <  q;  T.  - 


The   system  (2.0-.1)  can  now  be  written  as 


\ 


'IJ 


'2J 


,  1  1  J  <  q 


'pj 


V 


y 


:2.4)     .Bj$j.,  +  (I-Lj-R^l^   -  Tjlj^,  «  Sj   ,   2  <  J  <  q-l  ; 


q  q'^q 


lere  we  have  Introduced  the  identity  matrix,  I,  which  is 
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always  assumed  to  be  of  the  same  order  as  the  square  matrices 
to  which  It  Is  added. 

Further  simplification  Is  obtained  by  Introducing  the 

# 
(p  X  q)  -  dimensional  column  vectors   (or  q-dlmenslonal  com- 
pound vectors) : 


(2.5) 


li 


\ 


^      > 
^ 


L  J 
and  the   (pq)x(pq) 

matrices) : 
(2.6) 


'W 


(|) 


21 


^pq 
-  order  matrices  (or  qxq 


r     ■> 

/             S. 

Sl 

^1 

^2 

«21 

• 

= 

• 

• 

• 

• 

• 

• 

• 

\ 

V 

V          ) 

^      J 

^ 


o 


,  R 


order  compound 


6 


^ 


x 


0 


0 

B, 


0 


V 


J 


>> 


0   T 


1     0 
0  T^ 


0  T 


N. 


J 


♦The  vector  \  of  (2.5)  determines  an  "ordering"  [1]  of  the 
unknowns,  ^^\»   tout  this  order  need  bear  no  relationship  to  the 
sequence  In  which  the  Iterative  computations  are  carried  out. 
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The  system  of  linear  equations  (2.0-.1),  or  (2.4)  now 
becomes 

(2.7)  M$  -  (I-L-R-B-T)$  -  S  . 

Similar  formulations  may  be  Introduced  for  more 
general  boxmdarles  and  In  higher  dimensions.  In  par- 
ticular, if  the  boundaries  are  composed  of  coordinate 
segments,  the  matrices  L.  and  R>  remain  sqxiare  but 
of  different  orders  while  the  matrices  B^  and  T. 
become  rectangular.  Another  pair  of  "neighbors," 
({)...,,  appear  in  (2.0)  with  each  unit  increase  of  the 
dimension  and,  correspondingly,  additional  pairs  of 
matrices  must  be  introduced  in  a  manner  similar  to  those 
of  (2.6). 
5.   General  Single  Sweep  Iterations. 

We  svumnarlze  here  some  terminology  and  known  results 
for  a  class  of  Iterative  methods  for  solving  (2.7);  this 
class  is  defined  as  follows:  Let  the  coefficient  matrix, 
N,  be  written  as 

(5.0)  M  =  N  -  P  , 

where   |n|  /  0.  We  call  this  a  "splitting"  of  the  coeffi- 
cient matrix  and  the  system  (2.7)  becomes 
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(5.1)  N$  «  P$  +  S   , 


with  the  formal  solution 


(5.2)      $  -  (N-P)"^S  -  (I-N'^P)"V^S 


The  Iterative  procedure  Is  defined,  starting  from  some 
arbitrary  guess,  ^   ,  at  the  solution  vector,  by  the 
rec\irslon 

(3.5)      N$(^^  -  Pf(^-^)  +  S   . 

Thus  In  one  sweep  through  the  mesh  a  new  Iterate  Is  obtained. 
Applying  (5.5)  recursively  yields  for  the  v-th  Iterate 

(5.4)  l(^)-[I+(N-^P)+(N-^P)2+...+(N-^P)^^}N-V(N-^P)^^^^   . 

Thus  [2]  ^^^^  — >  J  as  V  — >  00  ,  for  arbitrary  "$}      ,   If  and 
only  If 

(5.5)  Llm  (N'-'-p)^  »  0   . 

V— >oo 

This  condition  Is  satisfied  provided  some  appropriate  norm 
[5]  of  (N~  P),  say  the  spectral  norm.  Is  <1. 


-  12  - 


This  result  la  more  frequently  obtained  by  Introducing 
the  sequence  of  error  vectors 


(3.6)         E^^^-l-]^^^^ 


which,  by  (5.1)  and  (3.5),  must  satisfy  the  homogeneous 
recursion 

(3.7)  NE^^^  -  PE^^'^^   ,   V  >  1   . 

The  eigenval\>es,  X,     of  N~  P  are  the   (pq)  roots  of 
the  characteristic  equation 

(3.8)  |XN  -  P|  »  0  . 

If  they  are  distinct*  there  then  exists  [4]  a  complete 
set  of  eigenvectors,  e.  ,  satisfying 

(3.9)  Xj^Nej^  .  Pej^   , 

which  span  the   (pq)  -  dimensional  vector  space.  Thus 
any  initial  error,  E^   ,  has  a  unique  expansion  in  these 
eigenvectors  of  the  form 


(3.10)         ,,,  JEl 


:(°)  =  H  a,.e. 


k«l 


♦  It  is  sufficient  here  to  assxnne  all  the  elementary  devisors 
[4]  of  (N"-^P)  are  simple. 
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By  (3.7)  and  (5.9)  the  above  yields,  for  the  v-th  error 
vector. 


Vk*k  • 
k»l 

In  order  that  ^^'^'   — >  ^  It  Is  necessary  and  sufficient 
that  E^^'  — >  0.  Thus  by  (5.11)  and  the  completeness  of 

the  eigenvectors,  convergence  Is  equivalent  for  arbitrary 

** 
Initial  error  to 

(5.12)         \   =  max  |>u  I  <  1   . 

k    *^ 

Let  us  require  that  the  most  slowly  decaying  component 

In  the  v-th  error,  (5. 11),  be  reduced  by  at  least  10~", 

where  m  >  0.  Then  we  must  have  y^   <  10~™  and  the  number 

of  Iterations  required  Is  bounded  by 

(5.15)  V  ^  -m/log  >v  =  m/R   . 

This  result  Is  valid  only  when  (5.12)  Is  satisfied  and  then 
R  a  -(log  A)  Is  called  the  rate  of  convergence.  This  quan- 
tity Is  useful  In  comparing  different  Iterative  methods  as 
the  number  of  Iterations  required  for  some  specified  con- 


♦*Por  a  special  Initial  error  In  which  the  amplitude  of  some 
particular  component  vanishes,  say  aj  =  0,  the  corresponding 
eigenvalue,  "K^,   Is  unrestricted.  However,  even  If  such  Initial 
distributions  could  be  determined,  computations  with  roundoff 
would  \andoubtedly  Introduce  this  component  at  an  early  stage 
of  the  Iterations  and  It  could  become  arbitrarily  large  If 
(5.12)  Is  violated. 
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vergence  criterion  varies  as  r"  . 


.-1, 


If  the  elementary  divisors  of  N  P  are  not  simple 
condition  (5.12)  still  suffices  for  convergence  but  the 
bound  (5.13)  does  not  apply.  To  obtain  such  a  bound  we 
assume  the  elementary  divisor  of  largest  order  correspond- 
Ing  to  X  »  |x  I   say.  Is  of  order  r  +  1.  Then  the  ex- 
pansion (5.10)  must  be  replaced  by  [4] 


(5.14) 


E 


(v) 


r+1 


r+l-k 

YL 

s»=0 


lYV'- 


s+k 


v-« 


k"k 


CH-2 


provided  v  >  (r+1)   (and  assuming  all  other  divisors  to 
be  simple).  The  largest  decay  factor  Is  now  given  by,  for 
V  >  (r+l)X  +  r  , 


v-r 


v-m 


and  to  reduce  the  error  by  at  least  10    reqxilres 
^v-r  I  M  <  10"™.   An  approximate  boiind  on  the  number  of  Iter- 
ations which  may  be  obtained  from  this  Inequality  Is 


(5.15) 


V  >  m/R  +  r  +  (r/2R)  log  [ (m/R  +  r)  m/R] 
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where  m  -  m-log  r!.  The  number  of  Iterations  required  for 
convergence  is  now  not  simply  proportional  to  R    but  this 
result  is  asymptotically  (for  m  — >  oo  )  equivalent  to 
(3.13).  For  practical  computations  the  discrepancy  may  be 
significant. 

A  large  rate  of  convergence  should  not  be  the  only 
criterion  in  evaluating  iterative  methods.  Rather  a  mea- 
sure of  the  time  required  (by  human  and/or  machine  effort) 
to  achieve  the  desired  accuracy  should  be  employed.  This 
time  is  essentially  proportional  to  the  total  number  of 
arithmetic  operations  (properly  weighted).  Thus  if  we 
let  N    be  the  operational  count  required  for  each 
iteration  (i.e.  to  solve  (3.3)),  the  total  time  is  pro- 
portional to  (assuming  (3.13)  to  be  valid), 

(3.16)         T=Nqp/R   . 

The  "best"  of  the  single  sweep  methods  is  that  for  which 
T  is  minimized.  Of  course  in  more  general  procedxares  a 
similar  criterion  should  be  employed. 

4.   Special  Single  Sweep  Iterations  and  Generalized 
ExtrapolatioiT 

With  the  above  notions  in  mind  we  consider  the 

special  class  of  splittings 
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(4.0)  N(t)  Z  V-'riL-YgR-TjB-T^T   ,  ?{y)   -  N(7)-M   ; 

where  the  real  numbers  y^     are  restricted  by  the  con- 
ditions that 

(4.1)  a)   |N(7)I  /  0,    b)   YiYgYjY^-  0  . 

This  defines  a  subset   P  of  the  5  dimensional  Euclidian 
space  of  all  points  y.     The  Iterations  are  defined  by 

(4.2)  N(y)  f^^^  -  P(7)$^''"'^^  +  S  . 

Thus  any  point  y€  f~  determines  an  Iterative  scheme 
given  by  (4.0)  and  (4.2);  we  shall  sometimes  refer  to 
this  as  the  scheme  y.     The  class  of  schemes,   f~.   Is 
Important  since  It  Includes  many  known  and  frequently 
used  schemes  while  the  data  arrangement  required  for  all 
of  these  schemes  Is  well  suited  for  automatic  computing 
machines.   Furthermore,  the  solution  of  the  system  (4.2) 
is  obtained  explicitly  when  7e  f""  in  one  sweep  over  the 
mesh  by  solving  either  two-term  or  three -term  recxirsions. 

In  particular  if  either  70*^1*^2  '^  ^  o**  '^o'^3'^4  ^  ^ 
three-term  recursions  are  Introduced  along  horizontal 
or  vertical  mesh  lines.   The  solution  of  these  recursions 
may  be  reduced,  by  the  well-known  factorization  of  a  Jacobl 
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matrix,  to  the  evaluation  of  two,  two-term  recursions 
along  the  appropriate  lines.  If  two  or  more  of  the 
7^  0  the  method  of  sweeping  the  mesh  to  solve  (4.2) 
In  one  sweep  Is  partially  determined  (I.e.  the  sweep 
must  start  at  a  particular  comer  or  with  some  line 
next  to  a  bo\mdary). 

The  operational  counts,  N  (7),  required  to  solve 
(4.2)  for  any  of  the  schemes  7€r~  vary  at  most  by  a 
factor  less  than  two.  The  minimum  niimber  of  operations 
needed  Is  4  multiplications  and  5  additions  at  each  mesh 
point  and  the  maximum  required  Is  7  multiplications  and 
6  additions  (neglecting  the  operations  done  only  once 
In  factoring  the  matrix  of  Jacobl  form).  The  special 
subclasses  of   f"  Introduced  In  Section  6  require  at 
most  6  multiplications  and  5  additions  at  each  point. 
Thus  In  the  remainder  of  the  paper  we  shall  ass\jme  the 
operational  coxmts  to  be  almost  equal  and  In  seeking  the 
"best"  scheme  we  shall  consider  only  the  eigenvalues. 

The  eigenvalues  of  a  scheme  such  as  (4.0-  .2)  are 
the  roots  X  of  the  characteristic  equation 

|XN(y)  -  P(7)l  -  0   ; 
or  explicitly,  of  the  eqviatlon 
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(^•5)        A  -  IgQl-g^L-ggR-g^B-g^Tl  -  0  ,  where 

Si  -  YiC^-l)  +  1  *   0  <  1  <  ^   • 

Each  root,  \.(t)*  is  ^^^s  si  function  of  the  scheme, 
•Y,  and  a  problem  naturally  suggested  Is  to  find  a 
y       €  n  such  that 

(4.4)      x(y»)  =  rain    {max\\{y)\)      . 

The  solution  of  this  problem  would  yield,  essentially, 
the  best  iterative  method  of  the  class  considered. 
The  usual  notion  of  extrapolation  (or  over-relaxation) 
of  the  scheme  y*     is  meaningless,  by  virtue  of  (4.4), 
since  no  improvement  on  rate  of  convergence  can  be  ob- 
tained.  In  fact  the  usual  extrapolation  techniques 
may  be  viewed  as  attempts  to  find  the  solution  of  (4.4) 
when  y     is  further  restricted  to  lie  on  some  particular 
curve  in   P  .  This  is  in  fact  shown  to  be  the  case  in 
section  6  and  the  curve,  in  the  usual  extrapolation  pro- 
cedures, is  a  straight  line.  The  problem  posed  in  (4.4) 
is  thus  a  generalization  of  extrapolation  in  which  the 
values  of  4  extrapolation  parameters  are  to  be  obtained 
(since,  by  (4,1b),  f     consists  of  four  4-dimensional 
subsets). 
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In  order  to  obtain  some  idea  of  the  possible  behavior 
of  "Kiy)     we  consider  schemes  y     on  the  "diagonal"  line 

These  points  do  not  lie  in  p  and,  obviously,  taking 
a  =  1  implies  direct  inversion  of  M.  However,  to  gain 
insight,  we  use  (4.5)  In  (4.5)  and  obtain  the  character- 
istic eqxiation: 

|[i(X-l)  +  1]M|  -  0   . 

Since,  by  assumption,   |M|  ^0,  there  is  only  one  eigen- 
value (with  multiplicity  pq)  and  it  is 

X  =  1-a 
Thus  along  the  diagonal  line  (4.5)  all  schemes  converge 
for  0  <  a  <  2  (i.e.  7j^  >  1/2)  and  diverge  otherwise 
(i.e.  y.   <  1/2).  Only  one  iteration  Is  required  for  a^l, 
as  then  X=0.  It  might  be  suspected  that  the  best  scheme 
in  P  is  obtained  by  taking  a  point  nearest  this  unit 
point  [a=l  in  (4.5)].  That  this  is  not  the  case  is  a 
simple  consequence  of  the  results  of  section  6.  The 
eigenvalues  of  the  above  example  become  arbitrarily 
large  in  absolute  value  as  a  — >  +  cd   and  the  point  y 
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approaches  the  origin;  the  roots  approach  unity  as 
a  — >  +  0  and  the  point  y     recedes  to  infinity.   It 
should  also  be  observed  that  these  eigenvalues  are  in- 
dependent of  the  matrix  M. 
5.  Some  Properties  of  the  Eigenvalues  X(7)« 

We  present  here  some  theorems  which  can  be  used  to 
compare  the  eigenvalues  of  various  schemes  y.     All  of 
these  results  are  simple  consequences  of  the 
Fundamental  Theorem;  Let  L,  R,  B  and  T  be  arbitrary 
matrices  of  the  form  (2.6),  (2.3).  Then  for  any  non- 
zero scalars  x  and  y, 

(5.0)         |M|  -  |I-L-R-B-T|  =  |l-xL-iR-yB-iT|   . 

—  —       X      y 

Proof:     Let  the  elements  of     M     be     M^„     where     1  <  (r, 
rs  —  \    ' 

s)  <  N  =  pq.  Then  each  term  in  the  formal  expansion  of 
|M|     is  given  by  [5] 

(5-1)  ±  \Mlf2M2) *'r,Tr(r) ^,rW      I 

where  t  is  one  of  the  NJ   permutations  of  the  first 
N  integers.   Let  each  point  (i,j)  of  the  original 
rectangular  mesh  (see  Fig,  1)  be  identified  with  a  unique 
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Integer  r  «  (j-l)p+l,  and  represent  any  permutation 
V    by  the  N  vectors  from  r  to  Tr(r)  on  this  mesh. 
By  the  definitions  (2.6),  (2.3)  and  (5.0),  M^,  /^v  ^  0 
only  If  ir(r)  =  r  or  the  point  corresponding  to  T(r) 
Is  an  adjacent  neighbor  of  the  point  r.  Thus  the  only 
permutations  which  lead  to  non-vanlshlng  terms  (5.1)  are 
those  whose  geometric  representation  Is  composed  entirely 
of  unit  vectors  In  the  (+1)-  and  (+J)  -  directions  and 
null  vectors.  However,  any  permutation  Is  a  product  of 
disjoint  cycles  [5]  and  the  representation  of  a  cycle 
Is  a  closed  path  on  the  mesh.  Thus  for  non-vanlshlng 
cycles  there  are  the  same  number  of  unit  vectors  In  the 
(+1)  -  direction  as  In  the  (-1)  -  direction  and  simil- 
arly for  the  (+J)  -  directions  (see  Fig.  1).   Since 
M   /  X   Is  an  element  of  L  If  ir(r)  =  r-1,  and  an 
element  of  R  If  T(r)  «  r+1,  there  are  as  many  factors 
from  L  as  from  R  In  each  non-vanlshlng  cycle.   A 
similar  result  Is  true  of  factors  from  B  and  T. 
Thus  each  non-vanlshlng  term  In  the  expansion  of  the 
right  hand  determinant  In  (5.0)  Is  Independent  of  x 
and  y  and  the  proof  Is  complete. 
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Figure  I. 


Geometric  representation  of  a  non-vanlshlng  cycle. 
In  the  permutation  of  which  this  cycle  Is  a  factor; 
ir(r)  «  r+p,  Tr(r+p)  =  r+2p,  ,  ir(r-l)  =  r. 


The  geometric  form  of  the  above  proof  was  suggested 
by  D.  Ludwlg.  Another  proof  has  been  given  by  B.  Friedman 
[8].  However,  the  geometric  proof  clearly  Indicates  that 
the  same  result  can  be  proved  for  a  matrix  which  represents 
the  difference  equations  on  any  connected  region  bounded 
by  coordinate  segments.   In  higher  dimensions  the  analog- 
ous theo.'em,  obtained  by  adding  two  more  matrices  to  M 
for  each  additional  dimension.  Is  easily  proved.   A  special 
case  of  the  above  theorem,  namely  that  obtained  by  setting 
X.  =  y  In  (5.0)  has  been  proven  by  D.  Young  [1].  Many 
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of  the  remaining  results  In  this  paper  hold  for  any 
matrix  M  which  can  be  written  In  the  form  (2.7)  such 
that  (5.0)  Is  satisfied  for  arbitrary  x  and  y. 

The  first,  almost  obvious,  consequence  of  the 
Fundamental  Theorem  Is  contained  In 

Theorem  I;  Let  X(7)  «  ><(7o''y2»'Y5»74)  ^e  any  eigenvalue 
of  the  scheme  y.     Then  for  all  7: 

(5.2)        >(7)  "  HyQ,y2*'^l''^J>''^k'^    "  ^(7o''Vi,'Y2»'Y4,73) 

Proof;  The  eigenvalues  X('y)  are  the  roots  of  the 
characteristic  eqviatlon  (4.3).  However,  by  the  Funda- 
mental Theorem  we  have 

A  =  IgQl-xg^L-  ^ggR-yg^B-  ^g^^Tl   , 


and  taking     x  -  ^gg/g^ /    y  -    ygi/gj^   the  above  yields* 

(5.3)  A  =.    IgQl-   /g^(L+R)-   yg^(B+T)|    . 

Thus  the  roots  of  A  =  0  are  Invar lent  under  the  Inter- 
changes 7^< — ^>  7^  ^"^  "^3  < — ^  ^k    '      Q'E.D. 

In  terms  of  Iterative  procedures  of  the  class  (4.2), 
Theorem  I  Indicates  that  various  ways  of  sweeping  the 


♦The  chosen  branches  of  the  sqxxare  roots  are  arbitrary  and 
hence  all  future  applications  of  (5.3)  could  be  written 
In  any  of  four  ways. 
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mesh  yield  identical  rates  of  convergence.   In  terms 
of  the  subset   f~  the  theorem  states  that  there  are 
certain  2-dimenslonal  planes  with  respect  to  which  the 
eigenvalues  are  symmetric.   Thus  the  volvime  of   f~ 
which  must  be  searched  for  a  best  scheme,  y* ,     is  re- 
duced by  a  factor  of  l/4.   Furthermore  it  is  of  interest 
to  note  that  along  the  2-dimensional  planes  y^   =  y^ 
and  y-,   =  yu     the  eigenvalues  ^(7)  must  have  relative 
maxima  or  minima  with  respect  to  the  appropriate  pair 
of  coordinates. 

A  somewhat  more  general  result  which  contains 
Theorem  I  as  a  special  case  is  contained  in 
Theorem  II;  Let  TCy)  be  some  particular  eigenvalue 
of  the  scheme  y.     Let  y     be  any  scheme  for  which  there 
exists  a  function  ?v(y)   such  that 

2 
(5.^)       /go]  &1&2       85^4   ,2  J  ^ 


^1^2   8384 

where  g^  and  g-   are  respectively  the  functions  of 
(T,  y)       and  (X,  7)  defined  in  (4.3).  Then  ^(7)   is 
an  eigenvalue  of  the  scheme  7. 

Proof;  Form  the  determinant  A  of  (4.3).  Then  as  in 
the  proof  of  Theorem  I  we  obtain  (5.3).   Using  (5.4) 
this  becomes 
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by  another  application  of  the  Fundamental  Theorem  to 
the  scheme  7.  However,  since  T{y)     is  an  eigenvalue 
of  7  we  have  A  =  0.  Thus  A  =  0  and  ^(7)  must  be 
an  eigenvalue  of  y.     Q.E.D. 

The  schemes  y     and  y,     and  the  eigenvalues  A(y) 
and  ^(7)  of  the  theorem  will  be  called  images  of  each 
other.   We  sometimes  refer  to  y     or  T{y)     as  the  refer- 
ence scheme  or  reference  eigenvalue  respectively.  It  is 
clear  from  {5 A)   that  the  relationship  of  being  images 
is  a  transitive  one.   Some  consequences  of  Theorem  II 
are  examined  in  the  remaining  sections. 

Another  theorem  which  is  qviite  useful  for  some 
special  difference  equations  is 

Theorem  III;  Let  the  matrices  (2.6)  be  such  that 
(L+R)  (B+T)  =  (B+T)  (L+r).  Then  as  is  well  known  these 
matrices  have  common  eigenvectors;   let  p.   and  m^  be 
the  eigenvalues  of  (L+R)  and  (B+T)  respectively, 
corresponding  to  the  common  eigenvector  e.  .  Then  each 
eigenvalue  of  any  scheme  7  is  a  root,  for  some  k,  of 
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the   (at  most  4th  degree)  eqioatlon* 

(5.5)       So  -  ygisj'pic  -  y8^^^k  =  °  • 

Proof;  Let  ^(7)  be  some  fixed  but  as  yet  unspecified 
eigenvalue  of  the  scheme  y.     Then  from  {h.J>),     A  =  0, 
and  as  In  the  proof  of  Theorem  I,  we  have  from  (5.3) 

(5.6)        IgQl  -  /g^(L+R)  -  yg^(B+T)|  =  0 

Thus  the  matrix  In  (5.6)  Is  singular  and  has  a  zero 
eigenvalue.  Applying  this  matrix  to  e.   we  see  that 

^k  -  s'o  -  ys^e?  Pk  -  JWk  ^\c 

Is  an  eigenvalue  belonging  to  the  eigenvector  e.  .  Using 
all  the  e.   we  obtain  all  of  the  eigenvalues  of  the  matrix 
In  (5.6).  However  at  least  one  i.    =  0.  Since  ^(7) 
was  any  eigenvalue  of  the  scheme  y     the  theorem  follows. 

This  theorem  Is  applied  generally  In  Section  7 
and  to  the  us\ial  Laplace  difference  equations  In  Section  8. 
6.  The  Complete  Image  Classes. 

The  pairs  of  Image  schemes  defined  by  Theorem  II  are  such 
that  only  one  eigenvalue  of  any  7  need  be  the  Image  of  one 
eigenvalue  of  the  reference  scheme,   7.   However,  of 

*As  pointed  out  In  the  proof  of  Thm.  I  there  are  four  such 
equations  which  could  be  used.  This  point  Is  clarified  in 
Section  7  where  more  of  the  consequences  of  the  hypothesis 
are  examined. 
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special  Interest  are  those  schemes  each  of  whose 
eigenvalues  Is  the  Image  of  a  corresponding  eigen- 
value of  some  particular  reference  scheme.  Such 
classes  of  schemes  are  called  complete  Image  classes 
and  we  proceed  to  obtain  three  of  them.  The  inter- 
sections of  these  classes  with  the  subset   I~  is 
grouped  into  5  subsets:   f"^,  Hg,  Hgi »  \~q     and 
Pc,,  of  which  Pg,  and  p^,  are  considered 
uninteresting. 

The  class  f..  In  order  that  any  two  schemes  y     and 
y     be  images,  (5.^)  must  be  satisfied  by  the  correspond- 
ing pair  of  image  eigenvalues.  However,  let  us  seek 
first  schemes  such  that 

(6.0)  s^S2  8384 

St  Sp       g-iSh 

is  an  identity  in  the  indeterminents  >  and  T.     This 
requires 

(6.1)  a)   71^2  =  7374      b)   7i72  =  73Y4 

Now  let  >  be  any  eigenvalue  of  a  scheme  y     which 
satisfies  (6.1b).  Then  if  7  is  any  scheme  satisfying 
(6.1a)  every  root  X  of 
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(6.2)     /gp,\  '    g^S^ 

is  an  eigenvalue  of  y     by  Theorem  II.  Since  (6.0) 
Is  satisfied  for  any  X  and  T,     each  equation  (6.2) 
obtained  for  a  different  eigenvalue  X  of  7  deter- 
mines at  least  one  eigenvalue  X  of  the  scheme  y. 
Thus  all  schemes  y     satisfying  (6.1a)  belong  to  the 
same  complete  Image  class,  say  class  A.  Since  the 
conditions  (6.1a)  and  (6.1b)  are  Identical  any  yek 
may  be  used  as  the  particular  reference  scheme.  We 
call  the  chosen  reference  scheme  7,  and  take  It 
to  be 

^A*   V  ^'   ^1-  Y2»  V  74-  0. 

This  Is  a  simultaneous  displacement  method  commonly 
known  as  Richardson's  method.  Thus  each  of  the  eigen- 
values of  any  class  A  scheme  Is  the  Image  of  one  of  the 
Richardson  eigenvalues. 

The  class  A  schemes  of  Interest  are  those  contained 
In  r~.  Prom  (4.1)  and  (6.1)  we  obtain  the  set  of  all 
such  schemes,   f"^,  and  they  are  listed  In  Table  A. 
The  real  parameters  a  and  3  are  to  be  restricted  such 
that  (4.1a)  Is  satisfied.  The  set   Pa  '^i-es   on  four 
2-dlmen8lonal  planes  In  p,  one  of  which  Is,  say: 
^2'=  Y4*  0  ,  7^«=  7^.   By  the  symmetry 
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•^0 

Yl 

-^2 

"^5 

•^4 

Ta 

1/a 

0 

1/p 

1/p 

0 

1/a 

1/P 

0 

1/p 

0 

1/a 

0 

1/p 

0 

1/p 

1/a 

1/p 

0 

0 

1/p 

^A 

1 

0 

0 

0 

0 

Table  A 

property  of  the  eigenvalues,  expressed  In  Theorem 
I,  we  need  examine  only  one  of  these  four  planes. 
However  using  any  76  P^  sii^cl  7^  In  (6.2)  yields 


(6.3) 


[^-(l-a)]2  =  t  ^A  t'^-(l-P)l   • 


This  equation  furnishes  the  mapping  of  the  eigen- 
values T.     of  7.   (Richardson)  onto  the  eigenvalues 
X  of  any  of  the  schemes  In  f"^. 

If  we  chose  a=«p»l  In  any  7e  f^  the  resulting 
scheme  Is  a  succusslve  displacement  method  commonly 
known  as  a  Llebmann  scheme.  In  the  present  connection, 
or  more  generally  as  the  Gauss-Seldel  method  (see 
Table  I).  The  mapping  (6.5)  yields 
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X  =  >k?,  or  X  =»  0  . 

Thus  some  T.  go  into  zero  and  others  into  their 
squares.  These  schemes  converge,  as  is  well  known 
[7],  if  the  Richardson  scheme  converges,  and  the 
rate  of  convergence  is  twice  as  large  (see  Section 

3). 

If  we  set  p  »  1  in  any  Yef";^  the  resulting 
scheme  is  a  successive  overrelaxation  [Ij  method 
where  a  is  the  overrelaxation  parameter  (see 
Table  I);  this  method  is  frequently  called  the 
extrapolated  Liebmann  scheme  [6].  The  schemes  7, 
in  this  case,  lie  along  four  lines  in  f". ,  for 
example:  Yg'^ii^^'  7l"73=l'  The  mapping  (6.5) 
becomes 

(6.5)        [X-l+a]2  =  a^y^X, 

which  has  been  studied  thoroughly  [1]. 

Taking  a  =  p  in  any  ye P^  yields  a  different 
successive  overrelaxation  method;  this  corresponds 
to  using  successive  displacements  (Liebmann)  over 
the  entire  mesh  and  then  extrapolating  (or  inter- 
polating) the  provisional  new  iterate  with  the  old 
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iterate.  The  mapping  (6.5)  becomes,  in  this  case, 
(6.6)       X  =  1  -  a(l-y?)  , 


which  is  easily  analyzed.  If  the  \  are  all  real 

A 

and  >.  <  1,  then  X  is  minimized  by  taking 


a  -  [1  -  ^&l  +^l)]'^ 


s/ 


where  X.  =  min|>. |.  This  method  is  the  analogue  of 
the  case  treated  at  the  end  of  Section  4;  it  is  called 
a  full-mesh  extrapolation. 

If  all  of  the  eigenvalues  T.  are 

A 

real  and  X.  <  1  the  general  mapping  (6.5)  can  be 

analyzed.  It  is  found,  in  this  case,  that  the  best 

of  the  class  A  schemes  is  extrapolated  Liebmann. 

However,  in  other  cases,  it  seems  possible  to  improve 

upon  the  Liebmann  extrapolations.  The  analysis  of 

this  mapping  has  been  done  by  K.  Cordis  and  will  be 

reported  in  a  future  paper. 

The  class  [""«.  Proceeding  as  in  the  previous  case 

we  now  require 

2 
/fo  I   ^  ^1^2 

1  Sq  I     ^182 
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to  be  an  identity  in  X  and  T.     This  is  equivalent  to 
a)  Yo-Yi-Tg  » 


to 
(6.7) 


b)  7o'7i»72 


Using  (4.1b)  these  relations  yield  the  indicated 
classes  Hg  and  Pg,  of  Table  B. 


•^0 

'>'l 

^2 

'^'3 

y^ 

Ib 

1/a 

1/a 

1/a 

1/& 

0 

1/a 

1/a 

1/a 

0 

1/3 

Fb' 

0 

0 

0 

1/a 

1/p 

Yb 

1 

1 

1 

0 

0 

Table  B 
The  canonical  reference  scheme,  Yg>  included  in  the 
table  is  a  simultaneous  line-displacement  method;  by 
its  analogy  with  7^  we  shall  call  it  a  line- 
Richardson  method.  The  classes  of  schemes   f^g  and 
r~g,  lie  on  three  2-dimensional  planes  in  f"  .  As 
before  if  we  take  X_  to  be  any  eigenvalue  of  y^ 
then  by  Thm.  II  any  root  X  of 


&r 


g. 
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Is  an  eigenvalue  of  7e  PU  (^^   f^B' ^  °^  ^^®  table. 
Of  course  the  appropriate  Y€ Pg  (or  ["ni )  ^s  to  be 
used  In  the  g..  Using  7€ Pg  and  7^  of  the  table 
this  yields 

(6.8)  [x.(i-a)]^  =^  ^I  [>^-(l-3)]  ; 

the  same  mapping  (6.3)  as  In  the   f".   schemes.  However, 
the  reference  schemes  7.  and  7^  are  not  Images  of 
each  other  and  thus  we  cannot.  In  general  compare  their 
eigenvalues.  In  Sections  7  and  8  special  cases  are 
considered  In  which  T.     and  7L     may  be  compared; 
It  Is  then  clear  that  some  schemes  In  f~~g  are  better 
than  the  best  scheme  In  f".  . 

All  of  the  simplifications  (6. 4-. 6)  are  shown  to 
apply  for  7€  f~n  ^y  choosing  the  same  special  a  and 
3  as  were  chosen  for  7e  ("a   (see  Table  I).  Thus 
we  may  call  y   =  (1,1,1,1,0)  "line-Llebmann" ,  and 
y   =  (l/a,  l/a,  l/a,  1,  O)  "extrapolated  llne-Llebmann" 
and  their  eigenvalues  have  exactly  the  same  relation- 
ships to  the  line -Richardson  eigenvalues  as  the 
eigenvalues  of  ordinary  Liebmann  and  extrapolated 
Liebmann  have  to  those  of  ordinary  Richardson. 

-  34  - 


The  remaining  class  of  schemes,  \~-Dtt     yield  the 
mapping 

(6.9)  [X-(l-a)][X-(l-3)]  «  aP^   . 

This  restilt  has  been  examined  and  Is  In  genez^l  fo\md 
to  be  Inferior  to  that  of  (6,8). 

The  class  f^Q   •  This  class  Is  determined  exactly  as 
the  previous  one  by  Interchanging  (7^,72^  with 
(7,,72^).  The  mappings  (6.8)  and  (6.9)  are  obtained 
for  the  corresponding  classes  l~'„     and  I~rt'     The 
canonical  reference  scheme  Is  y„   =  (1,0,0,1,1), 
a  line- Richardson  method,  which  Is  not  an  Image  of 
7^  or  Yg.   In  the  special  cases  of  the  next  two 
sections  It  Is  possible  to  compare  Xg  and  X-  and 
thus  to  determine  which  of  the  classes   pi  or   f". 
contains  the  best  scheme. 

The  essential  difference  between   [""„  and  f^ 
schemes  Is  the  direction  of  the  line  along  which  three- 
term  recursions  must  be  solved.   In  all   Pg  schemes 
they  are  parallel  to  the  direction  of  Increasing  1 
and  In  \~^     parallel  to  the  direction  In  which  J 
Increases. 
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7.  Reference  Eigenvalues  for  a  Special  Class  of 
Eqxiatlons. 

We  consider  here  those  systems  of  the  form  (2. 2-. 7) 
for  which  (L+R)  and  (B+T)  commute.  Then  the  Fund.  Thm. 
and  Thm.  Ill  are  valid.  However,  before  applying  these 
results  we  shall  examine  some  other  very  Important 
consequences  of  commutatlvlty  for  matrices  of  the 
Indicated  forms. 

Using  the  forms  (2.6)  In 

(L+R) (B+T)  =  (B+T) (L+R) 
we  obtain  the  conditions 

(Lj+Rj)Bj  =  Bj(Lj.^+Rj.^)   ,   2<  J  <  q   , 

(^j-^^j^^j  =  Tj(Lj+i-^«j+i)  '    1  <  J  <  q-1  ; 

These  conditions  are  necessary  and  sufficient  for 

* 
commutatlvlty.  With  no  loss  in  generality  we  may 

assxime      |B.|   f^  0     and      |T^|   ?^  0     in  which  case  the  above 

conditions  become 


(7.0)   (Lj+Rj)  = 


Bj(Lj.i+Rj.i)B-^  ,   2<  J  <q 


Tj(Lj^l+Rj^l)T-^  ,   1  <  J  <  q-1   . 


Since  the  (L.+R^)  are  all  similar  they  have  the  same 

♦Only  the  boundary  conditions  may  introduce  singular  BT 
and  Tj   for  elliptic  difference  equations.  However,  "^ 
these  cases  may  be  eliminated,  depending  on  the  boundary 
conditions,  by  either  not  including  the  boundary  values 
as  unknowns  (as  in  Section  8)  or  by  requiring  the 
difference  analogue  of  the  equations  to  hold  at  boundary 
points. 
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eigenvalues  and  elementary  divisors;  and  as  the  order  of 
these  matrices  is  pxp  there  are  at  most  p  distinct 
eigenvalues  p. ,  1  <  i  <  p.   Prom  the  form  (2.6)  of 
(L+R)  we  see  that  its  eigenvalues  are  Just  the  p., 
each  of  whose  multiplicity  has  a  factor  q  (i.e.  if 
p.   has  multiplicity  r.   in  (L^+R.)  then  it  has 
multiplicity  r^^q  in  (L+R)).   If  the  p.   are  all 
distinct  (or  belong  to  simple  elementary  divisors)  it 
can  be  shown  from  (7.0)  that  for  some  scalars  a^  , 

(7.1)  Bj^^  .  a^Tj^    ,   1  <  J  <  q-1   . 

By  changing  the  ordering  (2.2),  (2.5)  from  the 
original  row-ordering  to  a  column-ordering  (which  is 
equivalent  to  a  similarity  transformation  of  L,R,B  and 
T)  we  find,  as  above,  that  the  eigenvalues  m,.  of 
(B+T)  are  the  eigenvalues  of  p  similar  qxq  matrices. 
As  eigenvalues  of  (B+T)  each  \x.     has  a  multiplicity 
which  is  a  multiple  of  p  and  there  are  at  most  q 
such  eigenvalues.   It  is  well  known  that  commuting 
matrices  have  common  eigenvectors.  Thus  in  the  above 
case  there  exist  common  eigenvectors  e. .  such  that 

(7.2)  (L+R)e^j  =  Pj^e^j   ,   (B+T)e^j  =  ^^J®iJ 
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for  each  p.  and  \i.     belonging  to  different  eleraentai^r 
divisors  of  (L+R)  and  (B+T)  respectively. 

To  detennlne  some  properties  of  the  p^^  and  \x^ 
we  consider  the  characteristic  equations  of  which  they 
are  the  roots.  Exactly  as  In  the  proof  of  the  Frind. 
Thm.  we  see  that 


pI-(Lj+Rp|  =  |pI-xLj  -  |Rj 


for  arbitrary  x  ;^  0.  Thus  taking  x  »  -1  we  find 
that  If  p  Is  an  eigenvalue  of  (L.+R.)  then  so  Is 
-p.  Furthermore  If  p  Is  odd  there  Is  at  least  one 
zero  eigenvalue  and  additional  zero  eigenvalues  roust 
occur  In  pairs.  Analogous  results  are  true  of 
the  eigenvalues  \i     of  (B+T). 

We  are  now  In  a  position  to  apply  Thm.  III. 
We  select  first  7^  =  (1,0,0,0,0)  and  using  (7-2)  we 
may  replace  the  pair  (p^^-^M-j^)  °^   ^^®  theorem  by  the 
pair  (p^,Hj)  to  obtain  from  (5.5) 

(7.5A)         TijA  =  Pi  +  ^^J 


*The  properties  of  the  eigenvectors  (7.2)  and  the  sign 
parity  of  the  eigenvalues  explains  the  apparent 
ambiguity  pointed  out  In  Thm.  Ill,  In  which  any  of 
four  seemingly  different  equations  could  have  been 
used.  They  are  not  different  but  correspond  to  different 
labeling  of  the  eigenvalues.  ' 
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Similarly  using  y^   =  (1,1,1,0,0)  we  get 

and,  finally,  with  y^,   -  (1,0,0,1,1)  Thm.  Ill  yields 

We  have  introduced  an  obvious  double  subscript  notation 

for  the  T   .  These  are  all  of  the  roots  of  the 

corresponding  schemes  since  (7.3)  holds  for  all  p.  and 

M-j  .  Using  other  schemes  y     in  (5.5)  would  yield  many 

other  eigenvalues  explicitly.  However,  any  of  the  class 

r~^>  r~-Q     or  f~"p  scheme  eigenvalues  are  obtained  by 

using  (7.3)  in  (6.2)  etc. 

As  a  further  condition  let  us  assume  that  the  p^^ 

and  |Xj  are  all  real.   (A  sufficient  condition  for  this 

would  be  that  (L+R)  and  (B+T)  are  symmetric.  Then 

Bj^^  =»  T.  and  by  (7-1)  the  T.  and  B^  must  be  constants 

times  the  identity.  Similar  results  would  hold  for  the 

transformed  (L+R)  block  matrices).  Then  by  the  parity 

of  the  eigenvalues  p^  and  p..  (i.e.  since 

max|p.  I  =  max  p^  )  we  obtain  from  (70A) 
i   ^     i   -^ 

(7.4A)        T^-p+ii 
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Furthermore,  if  p  <  1  and  \i  <  1     we  get  from  (7.5B,C): 


(T-'lB)        ^  =  ^ 


* 


P 


Thus  we  may  easily  compare  these  principal  eigenvalues 
to  obtain 

a)  ^  +  ix  5  1  r==>  \  <:  {%,  \)     ; 


(7.5) 


b)  1^-1/2 1   $   |ir-l/2|   =>  Xb$  \ 


Since  the  mappings  of  the  eigenvalues  of  the  schemes  in 
|~  ,  f"  and  p^  from  the  corresponding  image  scheme 
eigenvalues  are  the  same,  the  best  scheme  will  be  foxuid 

A 

in  that  class  for  which  X  is  smallest.  The  relations 
(7.5)  may  thus  be  used  to  determine  the  best  class.  It 
should  be  noted,  from  (7.5a),  that  in  the  present  case 

A 

if  the  Richardson  method  converges  (X^  <  1)  then  the 
line-Richardson  methods  converge  faster,  and  if  Richardson 
does  not  converge  neither  do  the  line -Richardson  methods. 

8.  Laplace's  Equation. 
As  an  example  we  consider  V  0  =  0  in  the  rectangle 
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0<x<l,  0<y<l;  with  (j)  =  (given  function)  on 
the  boundary.  Using  the  mesh 


(8.0) 


and  the  difference  approximations 


<^x    Ax 

we  obtain  the  difference  equations 


(l<i<p-l 

(8.1)  0ij=e^[(|>i_i,j+1)i+i,j]+©y[<l>i,j.i+(^i,j+l],^  ^^^_^ 

Here  <})<-  i»  ^r.  ^»  ^i  n  ^^^     ^i  «  ^^®  ^^^  given 
boundary  values  and 

(8.2)  ©  = |2[i^  ,  9  = |xL^  ,  (e^   1). 

^   2(Ax%Ay^)      y   2(Ax2+Ay^)     x  y  2 

Equations  (8.1)  may  be  written  in  the  notation  of  Section 

2  where  i. .  =  r .  .  «  0  ,  b. .  =  t . .  =  9   and  the 
ij    ij    X*  ij    ij    y 

inhomogeneous  terms  come  from  the  boundary  values;  the 
system  is  of  order  [(p-l)(q-l)]  .  The  matrices  (L+R) 
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and  (B+T)  then  commute  and  are  symmetric,  and  as  proved 
in  Section  7,  their  eigenvalues  are  the  same  as  those 
of  respectively. 


r 


A 


9, 


0  1 

1  0  *•. 


•     • 


•   *    1 


•     • 


10 
""   (p-l)x(p-l) 


0  1 

1  0  **. 


A 


and   0. 


V 


1 

•   • 

1   0 

J 
(q-l)x(q-l) 


The  eigenvalues  of  these  matrices  are  easily  fo\md 


to  be 


Pi  '  ^®x  °°°  ^^  "^^^^    '     1  <  ^  <  P-^   * 

(8.5) 

M..   »  2©y  COS  (J  ir/q)  ,   1  <  J  <  q-1 

The  parity  properties  deduced  In  Section  7  are  seen  to 

be  satisfied  as  p^  ^   "Pp-i  ^"^  ^^   "  "^q-J*  "^^"^  ^^® 
principal  eigenvalues  are  obtained  for  l  =  j  =  1,  and 

we  have 

,2 


'p  =  2e„  cos  tt/p  ^  20  -  0  (tt/p) 


(8.4) 


/\ 


H  =  20y  cos  ir/q  ;ii  20y  -  ey(Tr/q)    , 
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where  the  approximate  values  hold  for  p  »  ir  and 
q  »  T. 

Let  us  consider  first  Iterative  solutions  of 
(8.1)  by  means  of  the  canonical  reference  schemes  of 
Section  6  (all  of  which  are  simultaneous  displacement 
methods).  The  principal  eigenvalues  of  these  schemes 
are  obtained,  since  Thm.  Ill  applies,  by  using  (8.4)  In 
(7«5);  then,  as  above,  retaining  only  terms  up  to 
second  order  In  l/p  and  l/q  we  have : 


[Richardson]; 

A  2   /©   e_\ 

(8.5)  y^^l-   Ig-I  ~f  "•■  ~i"  *  [Horizontal  Line -Richardson]; 
'^^'x  \^  p    q  / 

A      _2  /©    e_\ 

^p  «  1-  Is-  -|  +  -|  ) »  [Vertical  Llne-Rlchardson], 
y  \P    q  / 


In  this  approximation  we  see  that  the  direction  of  sweep, 

A 

for  minimum  "K     of  the  line-methods.  Is  determined  only 
by  the  mesh  ratio  and  does  not  depend  upon  the  number  of 
mesh  points  in  the  x  or  y-directlons  (since  both  p 
and  q  were  assumed  large).   If  the  nvunber  of  points  is 
"small"  the  analogous  criterion  is  obtained  by  using  the 
exact  eigenvalues  (8.4)  in  (7.5).  The  rates  of  convergence 
(Section  5)  of  the  above  schemes  are  easily  compared  by 
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recalling  that  -  log  (l-e)s;€  for  small  €.  If  the 

mesh  is  square  o  »  Q  »  l/4  and  the  line  methods 

A    y 

converge  twice  as  fast  as  the  ordinary  Richardson 

(to  second  order  In  l/p).  If  the  mesh  Is  rectangular 

then  ©  <  1/4  or  e_  <  1/4  and  by  sweeping  In  the 
^         y 

proper  direction  further  improvement  is  obtained  (the 
implicit  equations  should  come  from  lines  parallel  to 
the  direction  of  largest  mesh  spacing). 

The  successive  displacement  methods  corresponding 
to  the  above  reference  schemes  may  be  taken  as  (see 
Table  I) 

(8.6)  y^   »  (1,1,0,1,0),  7b  =  (1,1,1,1,0),  7g  »  (1,1,0,1,1) 

The  eigenvalues  of  these  schemes  are  related  to  the 
corresponding  reference  eigenvalues  (since  they  are 
complete  image  schemes)  by  (6.3)  with  a  »  3  =  1.  We  note 
as  in  Section  6  that  X  =  0  may  become  an  eigenvalue  of 

the  schemes  (8.6)  independently  of  the  values  of  T.     The 

—2 
remaining  roots  become  X  =■  X  ,  and  corresponding  to 

(8.5)  we  get,  up  to  second  order  in  l/p  and  l/q, 
\^  ^  l-2ir^    (  "^2"^  ^  )  *  [Liebmann]; 

2  /©    0 

(8.7)  Xg  <i^  1-  J-  (  -|  +  -J  )  >  [Horizontal  Line -Liebmann]; 


X(,  ^  1-  —■  (  -^  +  -^  )  »  [Vertical  Line-Liebmann]. 
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Of  coxirse  the  rates  of  convergence  are  now  rigorously 
twice  those  of  the  reference  scheme  rates.  The  best 
sweep  direction  Is  determined  as  In  the  previous  case, 
and  again  the  best  line  methods  converges  at  least 
twice  as  fast  as  the  ordinary  Llebmann  (to  second  order 
In  l/p  and  l/q). 

The  successive  overrelaxatlon  schemes  which  are 
usually  applied  to  the  above  are 

(8.8) 
YA-(l/a,l,0,l,0),7B-(l/a,l/a,l/o,l,0),'yc-(l/a,l,0,l/o,l/a), 

and  now  (6.5)  holds  In  each  case  with  p  «  1.   As  Is  well 
known  [1]  for  7.   above  (which  Is  extrapolated  Llebmann) 

A 

and  hence  for  all  these  cases,  "K     will  be  a  minimum 
when  a  Is  chosen  as  the  smaller  root  of 

(8.9A)  T^a^   -  4a  +  4  «  0  . 

[I.e.      Solve   (6.5)   for     "K     and  set  the  discriminant  eq\ial 

_    ^  A 

to  zero  when  X  =  X. ]  Then  we  have  all   |x|  »  X  where 

(8.9B)  X  »  a  -  1 

Using  the  reference  eigenvalues  (8.5)  In  the  above  we 
obtain  for  the  schemes  (8.8),  up  to  second  order  In 
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l/p  and  l/q, 
(8.10)       wg 

Xg»l-  ^2/9 J''^^t/-|+  -^  j  -^|-(-T"*^  1 ,  [ Ext . Horlz .  Llne-Lleb .  ] , 

X^=l-  ^2/9  V^TTJ -^4-|j  +|-j^-^  j.  [Ext.  Vert  .Llne-Lleb.  ] . 


There  Is  an  order  of  magnitude  Improvement  In  the  rates 
of  convergence  of  the  extrapolated  schemes  over  the 

previous  schemes.  That  Is,  since  (6^/p  +  ©y/Q  ) 

2  2   2   2  —1 
■  Ax  Ay  (Ax  +Ay  )   ,  the  present  rates  are  ^(Ax)  or 

^(Ay)  while  from  (8.5)  and  (8.7)  the  rates  are  ^(Ax^) 

or  (9^ Ay  ) .  The  line  schemes  now  Improve  the  convergence 

1/2 
rate  by  a  factor  2  '   for  a  square  mesh  or  larger  for 

rectangular  meshes  swept  In  the  proper  direction. 
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